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Abstract 



o 

L«. ' The field is expanded in a wavelet series and the wavelet coefficients are varied 

t::^ . in a simulation of the 2D (p^ field theory. The drastically reduced autocorrelations 

^> I result in a substantial decrease of computing requirements, compared to those in local 

Metropolis simulations. The improvement is shown to be the result of an additional 

freedom in the choice of the allowed range of change at the Metropolis update of wavelet 

Oh! components, namely the range can be optimized independently for all wavelet sizes. 
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1 Introduction 

In the present work we introduce a new approach to simulations of field theories, based 
on wavelet expansions. [^ [0] ||^ Wavelets form complete sets of localized, orthonormal 
states. Elements of a set of wavelets differ both in their locations and in their scales. In 
fact, complete sets contain wavelets at all scales, starting from an elementary scale up to the 
scale of the system. Wavelets can be labeled by two integers, n and fc, characterizing their 
scale and their location. They can be formed from a function ip^x) as 

^fc,„(x)=2-'^/V(2"x + A;). (1) 

Wavelets have been applied succcessfuUy to linear problems, like signal analysis. In the 
continuum, wavelets form a complete orthonormal system, allowing the expansion of L^ 
integrable functions, in a manner analogous to Fourier expansion. The orthogonal functions 
of the Fourier expansion are not localized, consequently they are not normalizable. Wavelets 
are localized, albeit on varying scales, from the elementary scale to the scale of the system. 
They are themselves L^ integrable function. As a rule, they lead to series expansions with 
convergence considerably faster then that of a Fourier expansion. 

The property of wavelets of analysing data at all scales makes them, in principle, at- 
tractive for the investigation of lattice problems. One can define wavelet representations, in 
which fields are expressed by their wavelet expansion coefficients, rather then by their values 
at given spacetime points. On finite lattices, the range of integers k and n in (|I]) is finite. It is 
not obvious that wavelet expansions, designed for investigating linear problems are useful for 
investigating a highly nonlinear problem such as a lattice field theory. The presence of large 
scale wavelets introduces nonlocality in the expression for interaction terms. The decrese 
in autocorrelation time is partially counterbalanced by the increased computational require- 
ments. It is then important to give estimates for the growth of nonlinearity in terms of the 
lattice size. The question is the balance: whether one gains more by decreasing autocorrela- 
tion times than one loses by having to deal with a number of nonlocal terms in the action. 
This question will be investigated extensively in the remainder of this paper. Wavelets have 
already been used in lattice problems: as a variational basis in the XY model 0], and as a 
tool to facilitate gauge fixing 0. 

Simulations in wavelet representation are somewhat similar to simulations using collective 
updating methods, such as Swendsen and Wang's algorithm |^ or the multigrid method [[7|. 
While the multigrid method is a so called smoothing filter, averaging a large number of 
points, the wavelet representation is, in a way, the opposite. Wavelets are designed to have 
their first few moments vanish. Collective updating methods have been designed to decrease 
the dynamical exponent z appearing in the expression of the autocorrelation time (r) of 
physical quantities 

r ^ c r, (2) 

where ^ is the spacetime correlation length, the inverse of the mass gap. In fact, some of the 
algorithms are able to reduce z to almost zero. Near second order phase transition points 
the correlation length is proportional to the linear size of the lattice, L. Thus, collective 
updating algorithms are particulary useful near critical points and for large lattices. In fact, 
for lattices of small or moderate size or at some distance away from the critical point, most 
collective updating algorithms are not much superior to local updating methods. 
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In the next section we will introduce the wavelet formalism and transform the 0^ field 
theory into wavelet representation. Then we will count the number of nonlocal terms, to 
get a handle on the increase of computational needs as a function of lattice size. In Section 
3. we put the wavelet method to a direct test, we simulate a 2 dimensional 0^ theory, both 
using the standard Metropolis algorithm and the wavelet method. The last section contains 
our conclusions. 



2 The Wavelet Method in Field Theories 

Wavelets form a orthonormal system of variables interpolating between coordinate and mo- 
mentum representations. In fact, wavelet representations (there are many different choices) 
unify the advantages of both. They are much more local and stable against small pertur- 
bations than the momentum representation, and at the same time they can describe long 
range correlations with a small number of terms, unlike the coordinate representation. 

The most widely used wavelet filters were discovered by Daubechies ||1|. Daubechies 
wavelets on discrete sets are defined as orthogonal matrix transforms. 2p of 2" spatial com- 
ponents, Xi,...,X2" are combined together in two different ways: with coefficients CiXi + 
C2X2 -l- ... + C2pX2p (smoothed combinations) and with different ordering and signs of the same 
coefficients —C2pXi + C2p-iX2 — ... + ciX2p (wavelets). The series is required to have p vanishing 
moments, i.e. J2k{~f^yc2p-k = 0, / = 0, 1, ...,p — 1. Similarly, one forms such linear combi- 
nation from subsequent components Xj, shifted by two lattice units an arbitrary number of 
times. Then requiring the orthogonality of the transformation makes the coefficients essen- 
tially unique. The wavelet transform consists oin—p subsequent orthogonal transformations 
of the above described type, applied to the smoothed combinations obtained in the previous 
orthogonal transformation. The coefficients Ci, ..., C2p are called the Daubechies wavelet filter 
coefficients. 

For p = 1 Co = Ci = I/V2, and the wavelets are called Haar- wavelets 0. The rth Haar 
wavelet of length 2^ on a one dimensional lattice is formed as 
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Xk,r = 2-'/' Yl [0(z + r2') - 0(z + 2'^-' + r2% (3) 



i=l 



where 0(j) denotes the field at the nth site. The range of k and r is k = 1,2, ...,n, and 
r = 0, ...,2"'"'^ — 1, respectively, where L = 2", is the size of the lattice. The orthogonal 
system is formed from the wavelets x{k,r), and one "smoothed" combination 

X = ^E0(^), (4) 

which is the normalized total magnetization. 

The structure of the wavelet coefficients for general p is similar to that of (^. Thus, the 
first 2"~^ wavelets are combinations of 2p neighboring data points shifted by multiples of 2. 
Then the next 2"~^ wavelets (each shifted from the next by 4 data points) are of length I2 = 
l + (22-l)(2p-l),etc. The2"-'=-i Mh type of wavelets are of length 4 = l + (2'=-l)(2p-l). 



They are shifted from each other by 2^ points. The last p combinations are not wavelets, 
but maximally smoothed combinations of the original data points. 

In an application to lattice field theory one has to use wavelets in more then one dimen- 
sion. In fact, wavelets can easily be defined in any number of dimensions. In each dimension 
the wavelets are labeled by their scale, ki, and by their location, r. Thus, a wavelet coefficient 
x{ki,ri; k2,r2; ...; kd^Td) corresponds to a wavelet box of length 2^^ in the ith. lattice dimen- 
sion, d is the number of dimensions. There are 2"'^~'^i~^2-...-fed wavelets labeled by a given 
scale vector k = {/ci, A;2, ..., kd}- The lattice field has a linear expansion in terms of the co- 
efficients x(^i;''"i; k2,r2] ...; kd,rd). We will call these coefficients the wavelet representation. 
It is on an equal footing with the momentum and coordinate representations. 

Wavelets are well suited to analyzing linear problems. Consequently a Gaussian model 
is just as simple in wavelet representation as in coordinate representation. Orthonormality 
implies that the form of the mass term is unchanged in wavelet representation 

^ = ^E%UOP (5) 

where for the sake of simplicity we have dropped the limits of the summations. The kinetic 
term is slightly more complicated, and its explicit form is not very enlightening. It has a 
moderate number of nonlocal terms. 

The propagator in wavelet space is quite analogous to the one in coordinate space. Thus, 
it has the same asympotic behavior and it also becomes singular at the critical point, -m^ = 0. 
Critical phenomena can be investigated in wavelet representation just as well as in coordi- 
nate representation. Alternatively, one can easily transform correlation functions back to 
coordinate representation, using the wavelet transformation. 

The question arises, however, how well the wavelet representation is suited to investigate 
nonlinear problems, like any nontrivial field theory. In the current paper we will investigate a 
typical example of nonlinear field theories, the 0*^ theory. We need to express the interaction 
term in terms of wavelets. Using the orthogonality of the wavelet transformation an explicit 
expression for the interaction term of the Lagrangian gJ24''^ is very complicated. It has 
many non-local terms, some of them correspond to 'long-range' interactions. The success of 
the wavelet representation depends on the average number of these long range interaction 
terms per wavelet. If this average number increases slowly with size we have a chance 
that simulations in wavelet representation will be more efficient then those in coordinate 
representation. 

First, we will calculate the average number of interaction terms per wavelet coefficient 
on a (i-dimensional lattice. It is sufficient to do the calculation in one dimension. If the 
number of interaction terms per wavelet is A^ in one dimension then it is N'^ in d dimensions. 
The interaction term is a quartic combination of the wavelet coefficients. For the sake of 
simplicity we will do the counting for the Haar-wavelet representation. Then it is easy to 
see that in every term the coefficient with the lowest value of j must be squared. Let that 
coefficient be x(j, m). Then the combinations of other wavelets multiplying xiJi^^Y can be 
found as follows: for each j there is only one wavelet overlaping x(j, -m). Then the bilinear 
multiplier has {n — j){n — j + l)/2 different types of terms. The number of wavelets with 



scale j is L/2K Then the total number of types of non-local interaction terms is 

^.. Y^L (n-j)in-j + l) 
Z^ of o 
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(6) 



In other words, the average number of non-local interaction terms per wavelet is 

(logL)2 



1^ 
2 2(log2)2 



(7) 



This number should be compared with 1, the average number of terms per site in the coor- 
dinate representation. In d dimension the appropriate factor is 



A^^ 



_(logL)2^ 
2d(iog2)2'^' 



In the next section, before discussing our simulations, we will perform a more precise 
analysis of computational needs. It will turn out that the number of computations per 
lattice site in a d dimensional model increases only as (logL)'' and not as (logL)^*^. The 
application of multiple hits at each wavelet coefficient in the Monte Carlo simulation further 
decreases the ratio of computational needs. 

3 Simulations 



The difference in the number of cycles needed between local and wavelet simulations is due 
to the calculation of the change of the action. The largest computational needs arise from 
the calculation of the interaction term. The action has the form 
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"^'0r + E(</'r-0r-e,)' + 



12^"^ 



(9) 



This expression is very complicated in terms of the wavelet coefficients, even for the simplest, 
Haar wavelets. For that reason it is not conducive to the calculations to express (|^) in 
terms of x(^i; ^ii ■■■ikdi ^d)- It is more efficient to store both the coordinate and the wavelet 
representation coefficients. This doubles storage needs, but accelerates and simplifies the 
code considerably. 

Suppose that during a sweep of the lattice we update the coefficient x = xi^i, ri, ...,kd, r^). 
Then the action has the following dependence on %: 
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+ terms independent of %, 



(10) 



The coefficients d(j)r/dx depend only on the relative position of the lattice points. They are 
given by products of coefficients q , for k = ki^k2, ■■■,kd. These coefficients can be easily 
tabulated with minimal storage requirements, even for general p. The coefficient of x^? and 
that of the quadratic part of the kinetic term depend on k only. For Haar wavelets the 
coefficient of the fourth order term is 2('^i+'^2+..+fed"dn)/2 

The action defines the probability distribution of x- It must be evaluated repeatedly 
during simulations. The calculation of the interaction part requires 8 floating point operations 
(5 for Haar wavelets) for every individual coordinate point, r. There are Ylhi — (2p — 
■]^^d2fci+...+fcd points inside the wavelet. Now, for a complete sweep one has to update all 
2nd-ki-...-ka coefficients of the same kind, giving the number of floating point operations 
as 8(2p — l)'^2"'°'. Finally, one has to update the coefficients for all values of ki,...,kd, 
bringing in an additional coefficient of n'^, where n = log L/ log 2, where L is the length 
of an edge of the lattice. Then the average number of floating point operations per lattice 
site is N^ = 8{2p — lYn'^. The calculation of the kinetic term requires further 4(2p — lYn"^ 
operations. For Haar wavelets the calculation of the kinetic term requires only 0{n'^^^) 
floating point operations. 

There is a further requirement in the updating algorithm. After every change of a variable 
X one has to recalculate the affected fields in coordinate representation. It is easy to get 
a similar estimate for the average number of floating point operations needed for such an 
updating. It just changes the coefficient of A'"^ from 8+4d to 10+ Ad. The ratio to the 
average number of floating point operations per site in a straightforward simulation in local 
simulation (coordinate representation) {Q+2d) is then 

t^ 10 + Ad, 



6 + 2d 



\n 



{2p-l)ni + 0{l/n)]. (11) 



where ti and t^ are the average times a sweep of the lattice takes in local and wavelet simu- 
lations, respectively. In particular, for Haar wavelets and two dimensions, Rd — (n^/2)[l + 
0{l/n)]. 

In fact, one can obtain a substantially smaller ratio as follows. Notice that the factor 
(logL)'^ appears in (pH]) due to the computation of the coefficients of x iii action (p^ . If 
one performs the Metropolis simulation in such a way that there are several hits at every 
wavelet x? then the same coefficients of (p!0D can be used and the extra calculation becomes 
not more costly then that of a hit in coordinate representation (local simulation). At the 
same time, having several hits at the same field component does not lead to any substantial 
savings in the local updating algorithm. Thus, in the limit of large number of hits the 
average sweeptimes for the two types of simulations converge. Of course, after a certain 
number of hits, the return in decreasing statistical errors diminishes. As usual, one has to 
find the optimal number of hits, so that the statistical errors were minimized using a fixed 
amount of computational time. Figure 1. shows the dependence of autocorrelation times 
on the number of hits. The autocorrelation times have been calculated using two different 
techniques. One, by direct measurement of the autocorrelation of the total magnetization, 
and the other by calculating the square of the relative error of the same physical quantity. 
They provide very similar information. When choosing the number of hits in the simulation, 
we had to keep in mind that the requirements on computing computing time slowly increase 
with the number of hits. The product of r^ with the average time a sweep, t^, should be 



optimized. We chose the optimum value of 5 in all of our simulations. For 5 hits we find 
the following values for t^/ti: 1.2, 1.3, 1.54, 1.82, on lattices of sizes L = 8, 16,32, and 64, 
respectively. 

The quality of an algorithm is determined by its improvement factor, which is defined as 

/=^, (12) 

where ti and r^ refer to real time autocorrelation lengths. Having determined (0) in the 
rest of this paper we will concentrate on the determination of autocorrelation times r^ and 

Tl- 

It is important to emphasize a substantial difference between a simulation using the stan- 
dard, local (coordinate representation) algorithm and a simulation in wavelet representation. 
A hidden parameter of all simulations is the size of the window. A, in which a new field value 
is picked during the Metropolis algorithm. In other words, if the value of the field component 
is 0, then we choose a random new value for 0, 0', in the interval (p — ^/2 < (p' < (p + A/2. 
While in a local simulation there is only one possible parameter A, in the wavelet repre- 
sentation there is the possibility, and indeed the necessity, of choosing different values Aj^, 
i,j = 1, ...,n + 1, for different sizes of wavelets. The effect of this freedom of choice of the 
window sizes on autocorrelation times will be discussed later. 

The evaluation of the improvement factor I of (|T^) requires the determination of auto- 
correlation times. We have performed extensive simulations on lattices of sizes L =8,16,32, 
and 64. In wavelet simulation, after a warmup period of O(IO^) sweeps we read out the 
relevant physical variables in the next O(IO^) sweeps. Using the local updating algorithm 
we had to have much longer runs, O(IO^) sweeps to be able to determine autocorrelation 
lengths. We compared the values of various physical quantities, such as (0^), (H), (|M|), 
{cv), and (x), where M is the magnetization, c„ is the specific heat, and x is the magnetic 
susceptibility. All the physical quantities agreed within statistical errors in the two simula- 
tions. The simulations were performed at fixed g = 1 and varying m^ between the limits of 
-0.31 and -0.25, a range which includes the critical point. 

The purpose of the simulations was to calculate autocorrelation times, obtained from the 
t dependence of autocorrelation functions 

r(t) = (x(t')x(t'-t)V, (13) 

where X is a physical quantity, t and t' label real time during the simulation, measured in 
units of sweeps. The average of ( |T3|) is over time t'. For moderate values of t, T(t) is expected 
to have an exponential dependence on t, T{t) ~ r(0) e"*/"^, where r is the autocorrelation 
time. In fact, the integrated autocorrelation function, 

1 * 1 

r^(t) = — - y r(t') ~ — . (i4) 

r(o) t'to 1 - e-^/" 

serves as a better measure of the autocorrelation time |^ |^ The integrated autocorrelation 
function goes into saturation at large t. The saturation value determines the autocorrelation 
time, r. 



We calculated autocorrelation times (measured in units of sweeps) for three physical 
quantities: Total magnetization, M, average value of 4>{x)^, and the average value of the 
action. The dependence of autocorrelation times on the mass parameter of the Lagrangian, 
near the critical point is shown in Figure 2. Typical errors are about 5-10%. The point of 
Fig. 2. represent data on lattice size L = 32. 

It is important to determine the dependence of autocorrelation times on the size of the 
lattice. Let us define an exponent C iii 8' manner, somewhat analogous to (0) 



T 



zL^. (15) 



At the critical point ( coincides with z. Away from the critical point ( is smaller then z 
because L > C,- The difference of z and ( measures the finite size effect on the correlation 
length. The exponent (, plotted in Fig. 3. in the critical region was calculated using all 
four lattice sizes we studied. Near the critical point, m'^ ~ —0.285, the exponent is almost 
2. In fact, if the smallest lattice, L = 8, is omitted from the fit the value exceeds 2. This 
shows that simulations in wavelet representation do not circumvent critical slowing down. In 
fact, a plot of exponent ( in local simulations is almost identical to Fig. 3. In other words, 
the possible gain in autcorrelation times (|T3p in local and wavelet simulations. Indeed, the 
ratio of autocorrelation times in local and wavelet representation simulations has a large and 
fairly constant value for all lattice sizes considered. Fig. 4. shows these ratios in the critical 
region of m^. A multiplier 5 corresponding to the number of hits per wavelet coefficient has 
been included in r^. The ratio has only a moderate variation with lattice size or m^. It is 



in the range of 10-15. Then the improvement factor, /, of (|12]) at the critical point varies 
between 7 on the largest lattice to a value over 10 on the smaller lattices. There is a slight 
tendency to an increase of / when one moves away from critical point toward the disordered 
phase. 

4 Discussion 

The general expectation of substantially decreased autocorrelation times in a simulation 
in wavelet representation was borne out by our simulations. One can demonstrate quite 
dramatically the effectiveness of our method on the free bosonic theory (^f = 0). Consider 
the autocorrelation of the total magnetization. As we mentioned earlier, the magnetization 
is proportional to one of the wavelet expansion coefficients, the 'non-wavelet' component. 
Let us denote this component by x- Since the kinetic term is independent of x, X completely 
decouples. The x dependence of the action is 

S = \m\\ (16) 

In other words, the action is a Gaussian function of x- Then the autocorrelation length for 
magnetization, can only depend on the dimensionless product rriA, where A is the range of 
change of x ^^ the simulation. The dynamical critical exponent is 2. This can be seen by 
the following heuristic argument. Suppose mA « 1. Then starting from an average value, 
X ~ 1/m, after n sweeps the value of x wanders away by an average value of 5x ~ a/^A 



(random walk). Complete decorrelation requires (^X ~ X ~ 1/m. Then we obtain the 
following estimate for the autocorrelation time: 

where A is the spatial correlation length on the lattice. ( p!7D shows not only that the critical 
exponent, z = 2, but also shows that the numerical value of the autocorrelation length is 
inversely proportional to A^. In other words, if one chooses A for the non- wavelet coefficient 
too small one can get very long autocorrelation times. In a wavelet simulation one can 
control A independently for different types of coefficients. Indeed, in our simulation we have 
chosen the window A for x in such a way that the acceptance was approximately 0.5. 

In a local, coordinate representation simulation one has no control over the allowed 
range of change of the total magnetization. In fact, it is expected to change by a very 
small amount in every sweep. Since X = Z) 0/vM, thus, 6x = Y^ 5(j)/yM^ roughly speaking 
5x ~ \fMb<\)l\fM = 5(f). We performed some simulations at (7 = as well (free bosonic 
theory), and in those simulations the optimal window A was independent of the mass, m. 
Since the autocorrelation can only depend on the product mA, at small m and constant A 
one can use the asymptotic behavior which should be the same as (p!7|). Indeed we obtained 
the following autocorrelation times for the total magnetization: 865 at m^ = .005, 310 at 
m^ = .01, 228 at m^ = .015, and 179 at m^ = .02. All of these have about 10% error. These 
autocorrelation times give a good fit to ([17|) . In other words, the effective value of A for 
the magnetization was indeed constant, as we assumed from the outset. Notice, however, 
what happens in a simulation in wavelet representation. The window A is adjusted for each 
type of wavelet separately. In particular, one adjusts A for the magnetization as well. If 
one requires constant acceptance, then Am and so the autocorrelation time is kept constant. 
The critical slowing down, at least for free theories, is completely eliminated. 

In the interacting case there is substantial mass renormalization. Still, for moderate 
values of the mA one expects a general dependence like (0) on A where one should use 
the renormalized mass. Unfortunately, in interacting theories in D = 2 there is another di- 
mensional quantity, g. The autocorrelation time also depends on the dimensionless quantity 
g/m?'. The dependence on this quantity seems to dominate the behavior of exponent C,. The 
constant c is, however, affected very strongly by the value of mA, leading to a dramatic 
decrease of autocorrelation times in simulations in wavelet representation. 

To check the above ideas, we ran a series of simulations in wavelet representation on an 
L = 32 lattice, at the same value of the physical parameters, set at the approximate location 
of the critical point, g = 1 and /x^ = —.285, but varying the window size for the wavelet 
component, that is proportional to the total magnetization. Fig. 5. shows autocorrelation 
times for magnetization as a function of the window size. A, for the appropriate wavelet 
component, x = Yl,<\>l\^ ■, where M is the total number of lattice sites. The curve has, for 
these values of m? and (7, a broad minimum around an optimal window size. 

The optimal value of window size for the change of in local simulations is about A = 3. 
Due to the orthogonality of the transformation from local variables, (/>, to wavelet compo- 
nents, X) ^ local simulation is in some sense equivalent to a wavelet simulation in which 
every window A is set equal to the value of A in the local simulation. For this reason we 
ran a simulation on a 32x32 lattice at m? = —.285, setting all windows equal to A = 3.0. 



We obtained an autocorrelation time of r^, = 226, comparable to r/ ~ 436. At the same 
point, if we optimize all the windows, such that all acceptances are between 0.4 and 0.5, then 
Tu, = 39. This result supports the idea that a major part of the decrease in autocorrelation 
time is not the result of relaxing the system at different scales simultaneously. It is rather 
the consequence of the extra freedom of setting the windows for every scale independently. 

In the present paper we have used Haar wavelets. It would be of some iterest to investigate 
simulations with other examples of the Daubechies wavelet series. Because the wavelet 
expansion coefficients would be all different, the calculation of the action would take a longer 
time. It is still possible that the improved convergence properties of the wavelet expansion 
would make the comparision with the local simulation even more favorable. 

Finally, we would like to point out that wavelet representations can be used for theories 
with more complicated order parameters as well. In such simulations the ratio of sweeptimes 
should be roughly the same as has been obtained in the current paper. As an example, for 
compact Abelian gauge theories wavelets should be formed from the components of the vector 
potential. The calculation of the change of the action requires the calculation of trigonometric 
function in both representations. The estimate we have given for computational requirements 
in Sec. 2. would still stand. Unfortunately, since the action is nonpolynomial, the ratio of 
simulation times may not be substantially improved by repeated hits at individual wavelet 
coefficients. 
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FIGURE CAPTIONS 

Fig. 1. 

Tijj for magnetization as a function of the number of hits per wavelet coefficient 
on a 32x32 lattice. The diamonds represent a direct measurement, while the 
triangles represent Tcfr calculated from relative errors. 

Fig. 2. 

T,uu as a function of m^ for the total magnetization (diamonds), for (0^) (triangles), 
and for the total energy (stars). 

Fig. 3. 

The exponent C as a function of rn^. 

Fig. 4. 

The ratio of autocorrelation times ti/t^j as a function of ni^ . Figures 4a, 4b, 4c, 
and 4d correspond to lattices of size L =8,16,32, qnd 64, respectively. 



Fig. 5. 



Autocorrelation time r^ as a function of A at rn^ = —.285 on a lattice of size 
L = 32. 



10 



